Consider the the number of pigs slaughtered in Victoria, available in the
aus_livestockdataset.
- Use the
ETS()function to estimate the equivalent model for simple exponential smoothing. Find the optimal values of \(\alpha\) and \(\ell_0\), and generate forecasts for the next four months.
## Series: Count
## Model: ETS(A,N,N)
## Smoothing parameters:
## alpha = 0.3221247
##
## Initial states:
## l[0]
## 100646.6
##
## sigma^2: 87480760
##
## AIC AICc BIC
## 13737.10 13737.14 13750.07
- Compute a 95% prediction interval for the first forecast using \(\hat{y} \pm 1.96s\) where \(s\) is the standard deviation of the residuals. Compare your interval with the interval produced by R.
auto_interval <- fc %>%
hilo() %>%
unpack_hilo("95%") %>%
select(4:7) %>%
slice(1)
auto_interval## # A tsibble: 1 x 5 [1M]
## .mean `80%` `95%_lower` `95%_upper` Month
## <dbl> <hilo> <dbl> <dbl> <mth>
## 1 95187. [83200.06, 107173.1]80 76855. 113518. 2019 sij
Let’s get manually those values:
sd_res <- fit %>%
augment() %>%
pull(.resid) %>%
sd()
auto_interval$.mean[1] + c(-1, 1) * (qnorm(0.975) * sd_res) ## [1] 76871.35 113501.77
Well, almost close …
Write your own function to implement simple exponential smoothing. The function should take arguments
y(the time series), alpha (the smoothing parameter \(\alpha\)) and level (the initial level \(\ell_0\)). It should return the forecast of the next observation in the series. Does it give the same forecast asETS()?
I’ve got an inspiration from this source. We are both getting the same values for this sub-exercise, but for the next one, this function yields better values.
exp_smoothing_manual <- function(y, arg_alpha, arg_ell_zero, bool_forecast_h1 = FALSE) {
if (bool_forecast_h1) {
total_len <- length(y) + 1
} else {
total_len <- length(y)
}
anti_alpha <- (1 - arg_alpha)
store_predictions <- rep(NA, total_len)
store_predictions[1] <- arg_ell_zero
for (i in seq_along(store_predictions)) {
if (i == 1) {
last_y <- store_predictions[i]
last_yhat <- store_predictions[i]
} else {
last_y <- y[i - 1]
last_yhat <- store_predictions[i - 1]
}
term_01 <- arg_alpha * last_y
term_02 <- anti_alpha * last_yhat
yhat <- term_01 + term_02
store_predictions[i] <- yhat
}
out <- yhat[length(yhat)]
return(out)
}
fit_model_pars <- coef(fit) %>%
select(term, estimate)
value_manual <- exp_smoothing_manual(
y = aus_pigs$Count,
arg_alpha = fit_model_pars$estimate[fit_model_pars$term == "alpha"],
arg_ell_zero = fit_model_pars$estimate[fit_model_pars$term == "l[0]"],
TRUE
)
value_auto <- fc$.mean[1]
value_manual == value_auto## [1] TRUE
Great!
Modify your function from the previous exercise to return the sum of squared errors rather than the forecast of the next observation. Then use the optim() function to find the optimal values of \(\alpha\) and \(\ell_0\). Do you get the same values as the
ETS()function?
Note: this source got value of 0.299 for alpha and 76379.265 for level.. But, the code below is heavily modified and doesn’t get the values of the source. Anyways, correct values are:
true_values <- coef(fit) %>%
select(term, estimate)
true_values## # A tibble: 2 × 2
## term estimate
## <chr> <dbl>
## 1 alpha 0.322
## 2 l[0] 100647.
Let’s try to get optimal values:
optimize_exp_smoothing <- function(pars = c(arg_alpha, arg_ell_zero), y) {
out_predicted <- rep(NA, length(y))
for (i in seq_along(out_predicted)) {
out_predicted[i] <- exp_smoothing_manual(
y = y[1:i],
arg_alpha = pars[1],
arg_ell_zero = pars[2]
)
}
squared_errors <- (out_predicted - y) ^ 2
out <- sum(squared_errors) %>% sqrt()
return(out)
}
optim_pigs <- optim(
par = c(0.5, aus_pigs$Count[1]),
y = aus_pigs$Count,
fn = optimize_exp_smoothing,
method = "Nelder-Mead"
)
true_values %>%
mutate(my_values = optim_pigs$par,
pct_diff = (my_values / estimate) - 1)## # A tibble: 2 × 4
## term estimate my_values pct_diff
## <chr> <dbl> <dbl> <dbl>
## 1 alpha 0.322 0.322 -0.0000323
## 2 l[0] 100647. 99223. -0.0141
Very small differences. For \(\alpha\), I am practically getting the correct value, while for \(\ell_{0}\) (starting value), I missed the mark by 1.41%.
Combine your previous two functions to produce a function that both finds the optimal values of \(\alpha\) and \(\ell_{0}\), and produces a forecast of the next observation in the series.
exp_smooth_combine <- function(time_series) {
optim_series <- optim(
par = c(0.5, time_series[1]),
y = time_series,
fn = optimize_exp_smoothing,
method = "Nelder-Mead"
)
out <- exp_smoothing_manual(
y = time_series,
arg_alpha = optim_series$par[1],
arg_ell_zero = optim_series$par[2],
TRUE
)
return(out)
}
var_forecast <- exp_smooth_combine(aus_pigs$Count)
var_forecast## [1] 95186.57
Is this same as forecasted value from fable?
## My Value: 95,186.57 | model value: 95,186.56